Critical Behavior in Graphene with Coulomb Interactions 
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We demonstrate that in the presence of Coulomb interactions, electrons in graphene behave like 
a critical system, supporting power law correlations with interaction-dependent exponents. An 
asymptotic analysis shows that the origin of this behavior lies in particle-hole scattering, for which 
the Coulomb interaction induces anomalously close approaches. With increasing interaction strength 
the relevant power law changes from real to complex, leading to an unusual instability characterized 
by a complex- valued susceptibility in the thermodynamic limit. Measurable quantities, as well as 
the connection to classical two dimensional systems, are discussed. 
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Graphene continues to fascinate physicists with its 
many unique properties 0, 0]. The low energy physics 
of electrons in graphene may be described by a Dirac 
Hamiltonian, analogous to the theory of massless neu- 
trinos, with the two components of the electron wave- 
functions representing amplitudes on the two sublattices 
that make up the graphene lattice. One class of prob- 
lems yet to be understood in a fuller description involves 
the effect of the Coulomb interaction. Its strength may 
be characterized by an effective fine structure constant 
P = e 2 /eHvp, where vf is the electron speed near the 
Dirac point and e is a dielectric constant due to a sub- 
strate. For suspended graphene, (3 is estimated to be of 
order 2. Naively, then, Coulomb interactions are rela- 
tively strong in graphene, so that one may expect to see 
its effects in a clean enough sample. This is the subject 
of our study. 

In the presence of Coulomb interactions, the Hamilto- 
nian of this system, remarkably, has no natural length 
scale: the 1/r interaction has precisely the same oper- 
ator dimension as the (Dirac) kinetic energy. This sug- 
gests that the system behaves as if it is at a critical point, 
even though no parameters need be tuned to attain this 
situation. In what follows, we demonstrate that a dra- 
matic effect of Coulomb interactions is that they induce 
power law correlations, a hallmark property of critical 
systems. The underlying cause of this originates in short 
distance physics - an anomalously large probability for 
close approaches of particle-hole pairs - but consequences 
are manifested at long distances because of the absence 
of a length scale in the Hamiltonian. The power law cor- 
relations in this system are reminiscent of the behavior 
of a variety of classical two-dimensional systems [3[ , and 
like those, when the coupling constant is sufficiently large 
we find indications of an unusual phase transition, char- 
acterized in the thermodynamic limit by a susceptibility 
that goes from real to complex rather than diverging. We 
speculate that resulting state may represent a precursor 
to the formation of a (gapped) exciton condensate 0, [B[ . 

The many-body physics we describe below is present in 



a simpler form for non-interacting Dirac electrons scat- 
tering from a Coulomb potential V(r) = —Ze 2 /r 0- 
H^ |. For small r, the wavefunctions vanish as ip m (r) ~ 
r V (m+i/2) 2 -z 2 /3 2 -i/2^ j n con trast, for impurity poten- 
tials which do not diverge so strongly at r = 0, the power 
law is fixed by the "centrifugal barrier" associated with 
a given angular momentum channel m, and does not de- 
pend on the potential itself. The fact that the exponent 
becomes a continuously varying function of /3 is unusual: 
the 1/r attraction allows an anomalous penetration of 
the centrifugal barrier. Moreover, for Z/3 > m, the expo- 
nent becomes complex, and one must introduce a short 
distance cutoff to obtain sensible wavefunctions. This 
"Coulomb implosion" phenomenon is accompanied by 
the appearance of a screening cloud p(r) ~ 1/r 2 around 
the impurity which is not present for smaller /?. 

General particle-hole channel propagators support 
analogous behavior, generating power law behavior in 
certain correlation functions. We show in particular that 
the sublattice-antisymmetric susceptibility does this, so 
that critical-like behavior is manifested in the response 
to potentials which are not symmetric for the two sublat- 
tices. Moreover we find a divergence when (3 exceeds /3 C , 
when the centrifugal barrier in particle-hole scattering is 
overcome and the power changes from real to complex, a 
many-body manifestation of Coulomb implosion. 

Analysis in Momentum Representation - To motivate 
our approach, we begin by analyzing the noninteract- 
ing problem in the presence of a Coulomb impurity as a 
scattering problem, using a standard momentum repre- 
sentation. The Hamiltonian is H a tp = [e— V(p)\ip, where 
Ha = hvpa ■ p is the kinetic energy for one of the two 
valleys, with a the Pauli matrices acting on the space of 
the two sublattices, p the momentum operator, and V(p) 
the Coulomb impurity potential. This is a low-energy 
effective Hamiltonian valid at distances large compared 
to the lattice scale; we ignore the small separation be- 
tween different sublattice points of the same lattice site. 
In clean noninteracting graphene there is negligible inter- 
valley scattering. This remains approximately true here 
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as well because the valleys are separated by a large mo- 
mentum, and the Coulomb interaction vanishes for large 
momenta. We only consider one spin species. The stan- 
dard (Lippman-Schwinger) equation for scattering states 
.13] 

in momentum space takes the form 



d 2 P ' 



(i) 

where is the unperturbed (matrix) Green's func- 

tion and ipw is an eigenstate for V = 0. Eq. ([l} is 
conveniently recast in terms of angular components, for 
which we define ip m (p) — J * ^e~" ne p^p(p), and decom- 
pose the Coulomb interaction in the form V(\p — p'l) = 
E„e- in ^-V)/ n (p'/p)/p, where 
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2vr [1 + x 2 - 2xcos(9)} 1 / 2 



(2) 



In terms of these quantities one arrives at a set of equa- 
tions coupling V'ltn an d V'am+i) wnere 1 , 2 are sub- 
lattice indices. To search for power law behavior in 
i)^ +s> {p) at large momentum p, we adopt the ansatz 

= Ca,m/p s for p — > oo. To lowest order in 1/p, 
non-vanishing solutions are supported when 



l-(Zf3) 2 I m ( S )I m +i(-s) = 0, 



(3) 



where I m {s) = / °° x 1 ^ 8 f- m (x)dx. Note that I m (s) di- 
verges for real s when s = 1 — |m| and s = 2 + \m\, 
and has a positive minimum value at s = 3/2. From 
these properties, one may see that Eq. © has solutions 
with real s when Zf3 is below a critical (Zf3) c ; above 
this s becomes complex. In particular, for m = 0, 
(Z[3) c = 1/2. This is Coulomb implosion expressed in 
momentum-space, which may be readily translated into 
a many-body context. 

Power Law Behavior for Interacting Dirac Electrons - 
The same analysis for a single particle-hole pair (the ana- 
log of the Cooper problem in superconductivity) reveals 
a short-distance power law changing from real to com- 
plex at a critical coupling [T3|. For the fully interacting 
many-body problem, such behavior should be contained 
in the four-point vertex function. Since the result is non- 
analytic in momentum, a nonperturbative approach is 
necessary. As we now demonstrate, the simplest relevant 
approximation scheme (a ladder sum (l5j ) shows a change 
in the exponent from real to complex in the particle-hole 
channel as (3 is increased. While a generic four point 
vertex includes the non-analytic behavior, the nonana- 
lyticity cancels in the density-density response By 
contrast, a potential that couples to the two sublattices 
asymmetrically (for example, by a substitutional impu- 
rity, or due to interaction with a substrate) , will include 
a response that is antisymmetric in the sublattice densi- 
ties. This antisymmetric response displays the power law 
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FIG. 1: (Color online) (a) Diagrammatical equation for the 
3-leg vertex T a a (k, q) with a z as the zeroth order vertex (the 
shaded cross in the figure), (b) Diagram for M(q). 



behavior. For the remainder of our discussion we focus 
on this antisymmetric response function. 

If the sublattice index is viewed as a pseudospin, the 
antisymmetric response can be viewed as a generalized 
magnetic susceptibility, defined as 



M{q) = --J dt([m z (-q,t),m z (q,0)}), (4) 
™>z(q) = o- z a/3 p a /3(q), p a p = ^ a \ + ^ a a k,^ ( 5 ) 



where A is the area of the sample, repeated indices are 
summed, and henceforward we will set K = 1. It can be 
shown that the reducible diagrams vanish in the compu- 
tation of M(q), so that we only need to consider irre- 
ducible bubble diagrams [Fig. l(b)| . 

The equation for the 3-leg antisymmetric vertex 
[Fig.|T(a)1 is 



d 2 q' 



(6) 



U{\q'\)K aMS {k - S)r#(fc - if, q), 



(2n) 2 

where U(\q\) = 2ire 2 /\q\, and 

K^&q) = i J ^GV>(p + q)G$(p). (7) 

Note that on the right-hand side of this expression, mo- 
menta are 3-vectors (go = 0), whereas elsewhere only the 
spatial components of the momenta remain. Defining 
X^(k,q) = K a ^5(k,q)T^(k,q), one finds 



(8) 



+K a f3~f8(k,q) 



d 2 q' 



U{W\) X ^(k-^q). 



This quantity is related to the susceptibility by 

d 2 k 



M(q) 



(2^ 



(9) 



Henceforth we focus on the long wavelength limit (small 
q), so for the moment we drop all terms of 0(q 2 ) and 
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higher. Using a circular moment expansion, one hnds 



~M(0) 
Aqq 



(k,q)=K^ B Jk,q)a 



-K 



(o) 

aafifl 



(k,q) 



aaj3 



(10) 



,k' 



Here we used the superscript (0) to denote the circular 
component m = 0, and the underlined indices are not 
summed over. In this equation we have introduced ex- 
plicit ultraviolet (A ~ 2n/a, a — lattice spacing) and 
infrared (fco ~ 27r/L, L = linear size of system) cutoffs. 

Defining x M (°)(fc, q) = a^xf^ (fc, q), in the limit 
q — » the solution to Eq. (fit)]) may be written in the 
form x M{O) (k,0) = ^F(^), where F obeys the inte- 
gral equation 



(in 



Note that F depends on the ratio fc/A, a reflection of the 
fact that the original Hamiltonian has no intrinsic length 
scale, so (in the limit fco — > 0) fc can enter only in this 
ratio. For fc/A <C 1, one easily confirms that Eq. (fTTj) is 
solved by a power law F(jh) ~ (A/fc) s , with s going from 
real to complex above some critical /3. 

To confirm this, we solved Eq. (|11[) numerically. For 
small (3, the solution is indeed a power law, provided fc ^> 
fco [see Fig. 2(a) inset]. For large enough /?, the solution 
is consistent with a power law of complex exponent, such 
that F becomes oscillatory with a power law envelope 
[Fig. [2(a)]. Moreover, M(q ->■ 0) = / -^X M {k,q ->■ 0) 
has a series of divergences [Fig. |2(b)| . The positions and 
weights of these poles depend on fc . We return to this 
important point momentarily. 

For small but nonzero q, it is interesting to compute 
the correction AM(q) = M (q) - M (0). The equation for 
the corresponding AF has a form very similar to Eq. (fTTj) , 
with only the "1" replaced by an inhomogeneous term, 
which is proportional to q 2 /k 2 for fc 3> q. The AM(q) 
resulting from this then vanishes with an exponent that 
varies with j3. The inset of Fig. [3] illustrates a typical re- 
sult for j3 not too large; the exponent as a function of /3 
is illustrated in Fig. [3l One physical consequence of this 
is that the difference in charge between sublattices for an 
impurity placed asymmetrically with respect to the sub- 
lattices will fall off with a /3-dependent power law at large 
distances, behavior which may be observable with a local 
scanning probe. We note that AM(q) has singularities 
at the same values of /3 as M (0), as should be expected 
from the form of Eq. (fTTj) . 

Figs. [2] and [3] are the central results of our work. They 
demonstrate that in the ladder approximation for the 
many-body problem of interacting electrons in undoped 
graphene: (i) For /? < j3 c generic particle-hole correlators 
decay with a power law at long distances, with an ex- 
ponent varying continuously with j3. The weak-coupling 
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FIG. 2: Solutions of Eq. JTTJ) with k /A = 10" 10 . (a) /3 = 0.5. 
Because the plotted is \F\, the oscillations appear as cusps. 
Note the amplitude of the oscillation scales roughly as 1/ \fk. 
Inset: F for /3 = 0.3. It is clearly a power law except for k 
close to fco ■ (b) The antisymmetric response M as a function 
of the interaction strength /3. 



many-body groundstate thus displays a basic property of 
a critical phase. The power law behavior is directly vis- 
ible in the sublattice-antisymmetric density correlator, 
(ii) For (3 > j3 c the exponent becomes complex, as in the 
noninteracting Coulomb implosion problem. In the in- 
teracting many-body case, the susceptibility M(q) of Eq. 
J5]) diverges for fc > 0. This strongly suggests a quantum 
phase transition to broken symmetry state with staggered 
charge order 0, 01- However, the presence of many such 
divergences as a function of /? suggests there are different 
ways to break the symmetry, (iii) The singularities van- 
ish in the thermodynamic limit, with the poles merging 
into a continuous function. The separation between them 
vanishes only logarithmically as fco — > 0, as we demon- 
strate below, resulting in a branch point at /3 C . We inter- 
pret this latter non-analytic behavior as the signal of a 
phase transition. Since it is a result of the merging poles, 
a natural interpretation is that the instability is into a 
state involving fluctuations among different realizations 
of a chiral order parameter which, if static, would pro- 
duce a gapped exciton phase [1, [f| . We speculate that 
with further increase in /3, one of these orderings could 
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FIG. 3: The exponent in AM(q) as a funciton of f3. Inset: 
AM(q) for P = 0.3. 



be favored over the others, resulting in a true condensed 
phase. 

Analytical Results for Model Kernel - A fuller under- 
standing of Eq. may be obtained by adopting a 
model kernel, 



fo(x) = 6(l-x) + -9(x-l). 

x 



(12) 



This has the same behavior as the real kernel at large and 
small x, and is simple enough to allow analytic solutions. 
We have verified numerically that the results for F and 
M are qualitatively very similar to those obtained with 
the correct /o- With this kernel, Eq. (jlip has general 
solutions of the form 



F(k) = A+k x + +A_k 



X- 



(13) 



-1±7 



and 7 = y/1 - 2/3. The 



with k = k/A, X± 
coefficients A± are determined by substituting Eq. (fT3|) 
back into the integral equation. This results in power 
law behavior for k 3> fco, with exponent A+, which goes 
from real to complex when f3 exceeds 1/2. Moreover, 
M(q —> 0) may be evaluated, yielding 



M(0) = 



A 



2k 



v F 1 + 7 - (3 + fc^(-l + 7 + (3) 
This has poles for f3 > 1/2 when 



V^2/3 — 1 In fco =2 arctan ^ - — h 2-7m, 



(14) 



(15) 



with integer n and < arctan (x) < it. Note that the 
distance between poles vanishes logarithmically as fco 
0, as discussed above. Furthermore, for j3 > 1/2, k^ 
becomes ill-defined unless an infinitesimal imaginary part 
is introduced in (3, so that (3 = 1/2 becomes a branch 
point for M(0). We interpret this as the signal of a phase 
transition in the thermodynamic limit, since M(0) need 
not be real and positive beyond this point. 



In this work we have examined the problem of interact- 
ing electrons in undoped graphene in the ladder approx- 
imation. For (3 < (3 C the ground state has /3-dependent 
power law correlations in the antisymmetric-sublattice 
density response. We note that in light of the the loga- 
rithmic growth ofvF 121 LL8| > one may expect logarithmic 
corrections to this power, which may be difficult to detect 
for finite size systems. For (3 > (3 C there is a nonstandard 
phase transition in which susceptibilities become diver- 
gent only for a finite size system. The continuously vary- 
ing exponent (in the ladder approximation) and the in- 
stability when it reaches a critical value are reminiscent 
of the behavior of two-dimensional classical XY mod- 
els. A central feature of our analysis is that, due to 
the absence of an intrinsic length scale, short distance 
power laws from anomalous penetration of a centrifu- 
gal barrier by the Coulomb interaction have an impact 
on the long-distance decay of correlators. Correspond- 
ing behavior occurs for non-interacting Dirac electrons 
near a Coulomb impurity. This analogy has been noted 
recently in a different wav(l9j]. 
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